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Abstract 

Large networks of spiking neurons show abrupt changes in their collective dy- 
namics resembling phase transitions studied in statistical physics. An example of 
this phenomenon is the transition from irregular, noise-driven dynamics to regu- 
lar, self-sustained behavior observed in networks of integrate-and-fire neurons as 
the interaction strength between the neurons increases. In this work we show how 
a network of spiking neurons is able to self-organize towards a critical state for 
which the range of possible inter-spike-intervals (dynamic range) is maximized. 
Self-organization occurs via synaptic dynamics that we analytically derive. The 
resulting plasticity rule is defined locally so that global homeostasis near the crit- 
ical state is achieved by local regulation of individual synapses. 

1 Introduction 

It is accepted that neural activity self -regulates to prevent neural circuits from becoming hyper- or 
hypoactive by means of homeostatic processes |[T4l . Closely related to this idea is the claim that 
optimal information processing in complex systems is attained at a critical point, near a transi- 
tion between an ordered and an unordered regime of dynamics J3j [TT] [9] . Recently, Kinouchi and 
Copelli [8 1 provided a realization of this claim, showing that sensitivity and dynamic range of a 
network are maximized at the critical point of a non-equilibrium phase transition. Their findings 
may explain how sensitivity over high dynamic ranges is achieved by living organisms. 

Self-Organized Criticality (SOC) 1 1 1 has been proposed as a mechanism for neural systems which 
evolve naturally to a critical state without any tuning of external parameters. In a critical state, typi- 
cal macroscopic quantities present structural or temporal scale-invariance. Experimental results |2| 
show the presence of neuronal avalanches of scale-free distributed sizes and durations, thus giv- 
ing evidence of SOC under suitable conditions. A possible regulation mechanism may be provided 
by synaptic plasticity, as proposed in iflOl , where synaptic depression is shown to cause the mean 
synaptic strengths to approach a critical value for a range of interaction parameters which grows 
with the system size. 

In this work we analytically derive a local synaptic rule that can drive and maintain a neural network 
near the critical state. According to the proposed rule, synapses are either strengthened or weakened 
whenever a post-synaptic neuron receives either more or less input from the population than the 
required to fire at its natural frequency. This simple principle is enough for the network to self- 
organize at a critical region where the dynamic range is maximized. We illustrate this using a model 
of non-leaky spiking neurons with delayed coupling for which a phase transition was analyzed in Q. 
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2 The model 



The model under consideration was introduced in [ 12] and can be considered as an extension of lfl5l 
5 1 . The state of a neuron i at time t is encoded by its activation level a; (t), which performs at discrete 
timesteps a random walk with positive drift towards an absorbing barrier L. This spontaneous 
evolution is modelled using a Bernoulli process with parameter p. When the threshold L is reached, 
the states of the other units j in the network are increased after one timestep by the synaptic efficacy 
Cji, di is reset to 1, and the unit i remains insensitive to incoming spikes during the following 
timestep. The evolution of a neuron i can be described by the following recursive rules: 



N 



a i(t) + J]] e ijHL(a,j(t)) + 1 with probability p 



Oi(f+l) = < 



N 



if di{t) < L 



a i{t) + 7 e ijHL(a>j(t)) with probability 1 — p 



N 

a .(t + l) = l+ tij H L{a 3 {t)) ifai(t)>L (1) 

where Hl{x) is the Heaviside step function: Hl(x) = 1 if x > L, and otherwise. 

Using the mean synaptic efficacy: (e) = J^f J2jj^i e ij/(N(N — 1)) we describe the degree of 
interaction between the units with the following characteristic parameter: 

V= w^wr (2) 

which indicates whether the spontaneous dynamics (r) > 1) or the message interchange mechanism 
(tj < 1) dominates the behavior of the system. As illustrated in the right raster-plot of Figure[T| at 
T] > 1 neurons fire irregularly as independent oscillators, whereas at rj = 1 (central raster-plot) they 
synchronize into several phase-locked clusters. The lower 77, the less clusters can be observed. For 
7] = 0.5 the network is fully synchronized (left raster-plot). 

In it is shown that the system undergoes a phase transition around the critical value r) = 1. 
The study provides upper (j max ) and lower bounds (T min ) for the mean inter-spike-interval (ISI) 
t of the ensemble and shows that the range of possible ISIs taking the average network behavior 
(At = T max -T m i n ) is maximized at 77 = 1. This is illustrated in FigureQ]and has been observed as 
well in |[8) for a similar neural model. 

The average of the mean ISI (r) is of order N x with exponent x — 1 for 77 > 1, x = 0.5 for 77 = 1, 
and x = for r/ < 1 as N — > 00, and can be approximated as shown in Q with [J 



3 Self-organization using synaptic plasticity 

We now introduce synaptic dynamics in the model. We first present the dissipated spontaneous 
evolution, a magnitude also maximized at rj = 1. The gradient of this magnitude turns to be simple 
analytically and leads to a plasticity rule that can be expressed using only local information encoded 
in the post-synaptic unit. 

3.1 The dissipated spontaneous evolution 

During one ISI, we distinguish between the spontaneous evolution carried out by a unit and the 
actual spontaneous evolution needed for a unit to reach the threshold L. The difference of both 
quantities can be regarded as a surplus of spontaneous evolution, which is dissipated during an ISI. 

'The equation was denoted (r)™,^ in (7). We slightly modified it using (e) and replacing 77 by Eq. (fJJ- 
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Figure 1: Number of possible ISIs according to the bound At = T max — T m i n derived in Q. 
For 77 > 1 the network presents sub-critical behavior and is dominated by the noise. For 77 < 1 it 
shows super-critical behavior. Criticality is produced at 77 = 1, which coincides to the onset of 
sustained activity. At this point, the network is also broken down in a maximal number of clusters 
of units which fire according to a periodic pattern. 

Figure^ shows an example trajectory of a neuron's state. First, we calculate the spontaneous evolu- 
tion of the given unit during one ISI, which it is just its number of stochastic state transitions during 
an ISI of length r (thick black lines in Figured). These state transitions occur with probability p 
at every timestep except from the timestep directly after spiking. Using the average ISI-length (r) 
over many spikes and all units we can calculate the average total spontaneous evolution: 

Etotal = ((t) - l) P . (4) 

Since the state of a given unit can exceed the threshold because of the received messages from the 
rest of the population (blue dashed lines in Figure |2^), a fraction of is actually not required to 
induce a spike in that unit, and therefore is dissipated. We can obtain this fraction by subtracting 
from (HJl the actual number of state transitions that was required to reach the threshold L. The latter 
quantity can be referred to as effective spontaneous evolution E e f f and is on average L — 1 minus 
(N — 1) (e), the mean evolution caused by the messages received from the rest of the units during an 
ISI. For 77 < 1, the activity is self-sustained and the messages from other units are enough to drive 
a unit above the threshold. In this case, all the spontaneous evolution is dissipated and E e ff = 0. 
Summarizing, we have that: 

If we subtract (|5]l from E to tai ©, we obtain the mean dissipated spontaneous evolution, which is 
visualized as red dimensioning in Figure^: 

E diss = E total -E eff = «r) - l)p - max{0, L — 1 — (N — l)(e)}. (6) 

Using (f3]l as an approximation of (r) we can get an analytic expression for Ediss- Figures^ and c 
show this analytic curve Edi SS in function of 77 together with the outcome of simulations. 

At 77 > 1 the units reach the threshold L mainly because of their spontaneous evolution. Hence, 
Etotal ~ E e ff and Edi SS ~ 0. The difference between E to tai and E e ff increases as 77 approaches 1 
because the message interchange progressively dominates the dynamics. At 7/ < 1, we have E e f f = 
0. In this scenario Edi SS — E to tai, is mainly determined by the ISI (r) and thus decays again for 
decreasing 77. The maximum can be found at 77 = 1. 

3.2 Synaptic dynamics 

After having presented our magnitude of interest we now derive a plasticity rule in the model. Our 
approach essentially assumes that updates of the individual synapses e,j are made in the direction of 
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Figure 2: (a) Example trajectory of the state of a neuron: the dissipated spontaneous evolution 
Ediss is the difference between the total spontaneous evolution E to tai (thick black lines) and the 
actual evolution required to reach the threshold E e f t (dark gray dimensioning) in one ISI. (b) Ediss 
is maximized at the critical point, (c) The three different evolutions involved in the analysis (pa- 
rameters for (b) and (c) are N = L = 1000 and p — 0.9. For the mean ISI we used r app of 
Eq. ©). 



the gradient of Edi SS - The analytical results are rather simple and allow a clear interpretation of the 
underlying mechanism governing the dynamics of the network under the proposed synaptic rule. 

We start approximating the terms N(e) and (N — 1) (e) by the sum of all pre-synaptic efficacies en-: 

N 

N(e) = (N- l)(e> + (e) « (N - l)(e> = e lk /N « ]T e ik . (7) 

i— 1 k^i k^i 

This can be done for large N and if we suppose that the distribution of is the same for all i. Edi SS 
is now defined in terms of each individual neuron i as: 



£ '''- = 1 ^ + V 1 *v 7 ~W 

-max{0,i-l-^e ife }. (8) 

An update of e,j occurs when a spike from the pre-synaptic unit j induces a spike in a post-synaptic 
unit i. Other schedulings are also possible. The results are robust as long as synaptic updates are 
produced at the spike-time of the post-synaptic neuron. 

Ae ij = K ^ = nf^L-^piL), (9) 

where the constant k scales the amount of change in the synapse. We can write the gradient as: 

~h ~ 



+ if(i-l- EM . eife )< 

, P J ---^indef if(L-l- E £ifc ) = 

^ ^(^t-^ + i) l-i if(i-i-E^^)>°- 

(10) 

For a plasticity rule to be biologically plausible it must be local, so only information encoded in the 
states of the pre-synaptic j and the post-synaptic i neurons must be considered to update . 
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Figure 3: Plasticity rule, (a) First derivative of the dissipated spontaneous evolution Edi as for 
« = 1, L = 1000 and c = 0.9. (b) The same rule for different values of c. 



We propagate J^k^i e ik to the state of the post-synaptic unit i by considering for every unit i, an ef- 
fective threshold L l which decreases deterministically every time an incoming pulse is received |6|. 
At the end of an ISI L l (L — 1 — Y^k^i €ik ) anc * encodes implicitly all pre-synaptic efficacies of i. 
Intuitively, L 1 indicates how the activity received from the population in the last ISI differs from the 
activity required to induce and spike in i. 

The only term involving non-local information in ( TTOb is the noise rate p. We replace it by a constant 
c and show later its limited influence on the synaptic rule. With these modifications we can write 
the derivative of E diss with respect to as a function of only local terms: 

dE dtss = -U-c + sgn(L') ^ 

deil 2yJ '{L 1 + 2c) 2 + 2c{L - L l ) 2 

Note that, although the derivation based on the surplus spontaneous evolution ( TTOb may involve 
information not locally accessible to the neuron, the derived rule (TUT i only requires a mechanism to 
keep track of the difference between the natural ISI and the actual one. 

We can understand the mechanism involved in a particular synaptic update by analyzing in detail 
Eq. (fTTT t. In the case of a negative effective threshold (L % < 0) unit i receives more input from 
the rest of the units than the required to spike, which translates into a weakening of the synapse. 
Conversely, if L 1 > some spontaneous evolution was required for the unit i to fire, Eq. (fTTb is 
positive and the synapse is strengthened. The intermediate case (L l — 0), corresponds to rj = 1 and 
no synaptic update is needed (nor is it defined). We will consider it thus for practical purposes. 

Figure^ shows Eq. ( fTTT i in bold lines together with dEl otal /deij (dashed line, corresponding to 
rj < 1) and dE l total / dtij + 1 (dashed-dotted, 77 > 1), for different values of the effective threshold 
L l of a given unit at the end on an ISI. E to tai indicates the amount of synaptic change and E e f f 
determines whether the synapse is strengthened or weakened. The largest updates occur in the 
transition from a positive to a negative U and tend to zero for larger absolute values of IS. Therefore, 
significant updates correspond to those synapses with post-synaptic neurons which during the last 
ISI have received a similar amount of activity from the whole network as the one required to fire. 

We remark the similarity between Figure [3p and the rule characterizing spike time dependent plas- 
ticity (STDP) HE). Although in STDP the change in the synaptic conductances is determined 
by the relative spike timing of the pre-synaptic neuron and the post-synaptic neuron and here it is 
determined by L l at the spiking time of the post-synaptic unit i, the largest changes in STDP occur 
also in an abrupt transition from strengthening to weakening corresponding to L % = in Figure^. 

Figure [3p illustrates the role of c in the plasticity rule. For small c, updates are only significant in a 
tiny range of IS values near zero. For higher values of c, the interval of relevant updates is widened. 
The shape of the rule, however, is preserved, and the role of c is just to scale the change in the 
synapse. For the rest of this manuscript, we will use c = 1. 
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Figure 4: Empirical results of convergence toward rj = 1 for three different initial states above (top 
four plots) and below (bottom four plots) the critical point. Horizontal axis denote number of ISIs 
of the same random unit during all the simulations. On the left, results using the constant k = 0.1. 
Larger panels shows the full trajectory until 10 3 timesteps after convergence. Smaller panels are a 
zoom of the first trajectory r/o = 1.1 (top) and rj = 0.87 (bottom). Right panels show the same type 
of results but using a smaller constant k = 0.01. 



3.3 Simulations 

In this section we show empirical results for the proposed plasticity rule. We focus our analysis on 
the time T conv required for the system to converge toward the critical point. In particular, we analyze 
how T conv depends on the starting initial configuration and on the constant k. 

For the experiments we use a network composed of N = 500 units with homogeneous L — 500 and 
p = 0.9. Synapses are initialized homogeneously and random initial states are chosen for all units 
in each trial. Every time a unit i fires, we update its afferent synapses e^ , for all j / i, which breaks 
the homogeneity in the interaction strengths. The network starts with a certain initial condition rjo 
and evolves according to its original discrete dynamics, Eq. (Q]l, together with plasticity rule (|9). 
To measure the time T conv necessary to reach a value close to r/ = 1 for the first time, we select a 
neuron i randomly and computer/ every time i fires. We assume convergence when r\ G (1 — v,l+v) 
for the first time. In these initial experiments, v is set to k/5 and k is either 0.1 or 0.01. 

We performed 50 random experiments for different initial configurations. In all cases, after 
an initial transient, the network settles close to r] = 1, presenting some fluctuations. These 
fluctuations did not grow even after 10 6 ISIs in all realizations. Figure [4] shows examples for 
t] Q e {0.58, 0.7, 0.87, 1.1, 1.3, 1.7}. 

We can see that for larger updates of the synapses (k — 0.1) the network converges faster. How- 
ever, fluctuations around the reached state, slightly above r] — 1, are approximately one order of 
magnitude bigger than for k = 0.01. We therefore can conclude that k determines the speed of 
convergence and the quality and stability of the dynamics at the critical state: high values of k cause 
fast convergence but turn the dynamics of the network less stable at the critical state. 

We study now how T conv depends on 770 in more detail. Given N, L, c and k, we can approximate 
the global change in 77 after one entire ISI of a random unit assuming that all neurons change its 
afferent synapses uniformly. This gives us a recursive definition for the sequence of 774 s generated 
by the synaptic plasticity rule: 



A(th) - k(N - 1) 



-LeffiVt) - c , sgn(?7 t - 1) 



! \j (LeffM + 2c) 2 + 2c(L - L eff ( Vt )) 
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Figure 5: Number of ISIs (a) and timesteps (b) required to reach the critical state in function of 
the initial configuration r/o- Rounded dots indicate empirical results as averages over 10 different 
realizations starting from the same t]q. Continuous curves correspond to Eq. ( fT2b . Parameter values 

are N = 500, L = 500,p = 0.9, c = 1, v = k/5. 



where L e //(»7t) = (i - 1) I 1 - -J and r) t +i = rjt + A(?yt). 
Then the number of ISIs and the number of timesteps can be obtained b>Q: 

Tconv = min({i : \<r) t - 1| < i/}), T co ™. s t e ps = X! T app^t)- ( 12 ) 

t=o 

Figure [5] shows empirical values of t cotlv and T conv _ s t e ps for several values of 770 together with the 
approximations 112) . Despite the inhomogeneous coupling strengths, the analytical approximations 
(continuous lines) of the experiments (circles) are quite accurate. Typically, for rj < 1 more spikes 
are required for convergence than for 770 > 1. However, the opposite occurs if we consider timesteps 
as time units. A hysteresis effect (described in Q) present in the system if 770 < 1, causes the 
system to be more resistant against synaptic changes, which increases the number of updates (spikes) 
necessary to achieve the same effect as for 770 > 1. Nevertheless, since the ISIs are much shorter for 
supercritical coupling the actual number of time steps is still lower than for subcritical coupling. 



4 Discussion 



Based on the amount of spontaneous evolution which is dissipated during an ISI, we have derived a 
local synaptic mechanism which causes a network of spiking neurons to self-organize near a critical 
state. Our motivation differs from those of similar studies, for instance fl8], where the average 
branching ratio a of the network is used to characterize criticality. Briefly, a is defined as the 
average number of excitations created in the next time step by a spike of a given neuron. 

The inverse of 77 plays the role of the branching ratio a in our model. If we initialize the units 
uniformly in [1, L], we have approximately one unit in every subinterval of length rje, and in conse- 
quence, the closest unit to the threshold spikes in I/77 cases if it receives a spike. For 77 > 1, a spike 
of a neuron rarely induces another neuron to spike, so a < 1. Conversely, for 77 < 1, the spike of a 
single neuron triggers more than one neuron to spike (a > 1). Only for 77 = 1 the spike of a neuron 
elicits the order of one spike (er = 1). Our study thus represents a realization of a local synaptic 
mechanism which induces global homeostasis towards an optimal branching factor. 

This idea is also related to the SOC rule proposed in [3], where a mechanism is defined for threshold 
gates (binary units) in terms of bit flip probabilities instead of spiking neurons. As in our model, 
criticality is achieved via synaptic scaling, where each neuron adjusts its synaptic input according to 
an effective threshold called margin. 

2 The value of T app (rit) has to be calculated using an (e) corresponding to rjt in Eq. ^3). 
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When the network is operating at the critical regime, the dynamics can be seen as balancing between 
a predictable pattern of activity and uncorrelated random behavior typically present in SOC. One 
would also expect to find macroscopic magnitudes distributed according to scale-free distributions. 
Preliminary results indicate that, if the stochastic evolution is reset to zero (p = 0) at the critical 
state, inducing an artificial spike on a randomly selected unit causes neuronal avalanches of sizes 
and lengths which span several orders of magnitude and follow heavy tailed distributions. These 
results are in concordance with what is usually found for SOC and will be published elsewhere. 

The spontaneous evolution can be interpreted for instance as activity from other brain areas not 
considered in the pool of the simulated units, or as stochastic sensory input. Our results indicate 
that the amount of this stochastic activity that is absorbed by the system is maximized at an optimal 
state, which in a sense minimizes the possible effect of fluctuations due to noise on the behavior of 
the system. 

The application of the synaptic rule for information processing is left for future research. We ad- 
vance, however, that external perturbations when the network is critical would cause a transient 
activity. During the transient, synapses could be modified according to some other form of learning 
to encode the proper values which drive the whole network to attain a characteristic synchronized 
pattern for the external stimuli presented. We conjecture that the hysteresis effect shown in the 
regime of r\ < 1 may be suitable for such purposes, since the network then is able to keep the same 
pattern of activity until the critical state is reached again. 
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